R 중급 통계분석 4일차 (2)

지도 시각화

1. maps를 이용한 시각화

maps 패키지

• 전 세계 나라의 경계선이 입력된 지도 자료이다.

• 나라의 주별 경계선까지 입력되어 있으나 지도가 1990년 자료를 기준으로 개발되어 업데이트 된 나라와 업데이트 안된 나라가 섞여있다.

• col과 fill을 통해 나라별 특징을 시각화 하기 좋다.

• 미국의 경우 state 까지 세분화 하여 시각화가 가능하다.

library(maps)
library(ggmap)
map("world")

map("world","italy")

map(‘world’)를 이용하면 전세계지도

만약 나라의 이름을 추가한다면 이에 해당하는 나라지도를 그릴 수 있다.

map_data() 함수를 이용하여 원하는 나라에 해당하는 자료를 가지고 올 수 있습니다만, 내용이 빈약합니다.

korea <- map_data("world","South Korea")
head(korea, 20)

##        long      lat group order      region subregion
## 1  126.3270 33.22363     1     1 South Korea  Cheju Do
## 2  126.2820 33.20151     1     2 South Korea  Cheju Do
## 3  126.2402 33.21484     1     3 South Korea  Cheju Do
## 4  126.2290 33.22524     1     4 South Korea  Cheju Do
## 5  126.1787 33.28257     1     5 South Korea  Cheju Do
## 6  126.1656 33.31201     1     6 South Korea  Cheju Do
## 7  126.1994 33.36806     1     7 South Korea  Cheju Do
## 8  126.3377 33.46040     1     8 South Korea  Cheju Do
## 9  126.6955 33.54932     1     9 South Korea  Cheju Do
## 10 126.7599 33.55322     1    10 South Korea  Cheju Do
## 11 126.9012 33.51514     1    11 South Korea  Cheju Do
## 12 126.9313 33.44385     1    12 South Korea  Cheju Do
## 13 126.9054 33.38237     1    13 South Korea  Cheju Do
## 14 126.8729 33.34116     1    14 South Korea  Cheju Do
## 15 126.7092 33.27168     1    15 South Korea  Cheju Do
## 16 126.5817 33.23833     1    16 South Korea  Cheju Do
## 17 126.3270 33.22363     1    17 South Korea  Cheju Do
## 19 126.7530 34.34399     2    19 South Korea        11
## 20 126.7699 34.29644     2    20 South Korea        11
## 21 126.6891 34.30542     2    21 South Korea        11

unique(korea$subregion)

##  [1] "Cheju Do"  "11"        "3"         "21"        "5"         "Namhae Do"
##  [7] "Koje Do"   "8"         "Ullung Do" "10"        NA

한국지도를 시각화 한 후 polygon 함수를 이용하면 제주도에 색을 넣을 수 있지만,

그림의 해상도가 낮아 시각화에 적절한 그림은 아닙니다.

map('world','South Korea')
jeju <- subset(korea, korea$subregion=="Cheju Do")
polygon(jeju$lon,jeju$lat, col="grey")

하지만, 미국의 지도는 maps 패키지만으로 시각화가 가능합니다.

map("state")
map.text("state", add=T)

2. raster를 이용한 시각화

raster 패키지는 인터넷을 통해 행정, 해발도, 기후 자료들을 다운로드 받아 이를 시각화 할 수 있습니다.

시각화를 위해선 raster 패키지의 설치가 필요합니다.

한국지도를 시각화 한다면 아래와 같이 코드를 입력하면 됩니다.

library(raster)
kor0 <- getData('GADM', country='KOR', level=0)
plot(kor0)

kor1 <- getData('GADM', country='KOR', level=1)
plot(kor1)

kor2 <- getData('GADM', country='KOR', level=2)
plot(kor2)

level 수준에 따라 0은 나라, 1은 시도, 2는 시군구 행정자료입니다.

getData 함수로 다운로드 받은 지도는 plot 함수를 이용하면 간단히 시각화 할 수 있습니다.

이제 raster 패키지에서 다운 받은 그림자료를 이용하여 코로나바이러스감염증-19의 시도별발생동향을 시각화 하겠습니다.

최근 시도별 발생동향

시도별 발생동향을 시각화 하기 위해선 아래의 단계가 필요합니다.

1. kor1 자료에 발생 현황 수 입력

2. 입력된 자료 범위 선정 후 색 부여

3. col을 이용한 지도 시각화

그럼 차례대로 지도를 그려보겠습니다.

class(kor1)

## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"

head(kor1@data)

##    GID_0      NAME_0   GID_1            NAME_1
## 1    KOR South Korea KOR.1_1             Busan
## 10   KOR South Korea KOR.2_1 Chungcheongbuk-do
## 11   KOR South Korea KOR.3_1 Chungcheongnam-do
## 12   KOR South Korea KOR.4_1             Daegu
## 13   KOR South Korea KOR.5_1           Daejeon
## 14   KOR South Korea KOR.6_1        Gangwon-do
##                                                                                                  VARNAME_1
## 1                                                          Pusan|Busan Gwang'yeogsi|Pusan-gwangyoksi|Fusan
## 10 Chungchongbuk-Do|Chungcheongbugdo|Ch'ungch'ong-bukto|Chusei Hoku-do|North Chungchong|Ch'ungch'ong-bukto
## 11                   Chungchongnam-Do|Ch'ungch'ong-namdo|Chusei Nan-do|South Chungchong|Ch'ungch'ong-namdo
## 12                                                        Taegu|Daegu Gwang'yeogsi|Taegu-gwangyoksi|Taikyu
## 13                                             Daejeon Gwang'yeogsi|Taej<U+014F>n-gwangy<U+014F>ksi|Taiden
## 14                                               Kang-Won-Do|Gang'weondo|Kangwon-do|Kogen-do|South Kangwon
##                  NL_NAME_1       TYPE_1          ENGTYPE_1 CC_1 HASC_1
## 1  부산광역시 | 釜山廣域市  Gwangyeoksi  Metropolitan City <NA>  KR.PU
## 10     충청북도 | 忠淸北道           Do           Province <NA>  KR.GB
## 11     충청남도 | 忠淸南道           Do           Province <NA>  KR.GN
## 12 대구광역시 | 大邱廣域市  Gwangyeoksi  Metropolitan City <NA>  KR.TG
## 13 대구광역시 | 大邱廣域市 Gwangyeoksi  Metropolitan City  <NA>  KR.TJ
## 14         강원도 | 江原道           Do           Province <NA>  KR.KW

class 함수를 실행시키면 kor1의 속성을 보실 수 있습니다.

공간자료로 구성되어 있으며, 공간자료는 @(slot) 개념을 이용하여 데이터를 관리하고 있습니다.

자료를 보면 한국이름 중 대전이 대구로 잘못 표기되어 있습니다.

따라서 NAME_1을 이용하여 자료를 처리하겠습니다.

이제 NAME_1의 순서대로 확인자수를 확인하여 입력해주세요.

kor1@data$NAME_1

##  [1] "Busan"             "Chungcheongbuk-do" "Chungcheongnam-do"
##  [4] "Daegu"             "Daejeon"           "Gangwon-do"       
##  [7] "Gwangju"           "Gyeonggi-do"       "Gyeongsangbuk-do" 
## [10] "Gyeongsangnam-do"  "Incheon"           "Jeju"             
## [13] "Jeollabuk-do"      "Jeollanam-do"      "Sejong"           
## [16] "Seoul"             "Ulsan"

confirmed <-c(26151,10928,16287,14350,8309,9706,11173,102984,14540,23941,20455,3900,11259,11454,2514,66067,8302)
kor1@data$confirmed

## NULL

이제 입력된 숫자를 cut 함수와 heat.colors를 이용하여 색을 넣겠습니다.

cut 함수는 임의의 구간으로 자료를 잘라 줍니다.

cut(confirmed, c(0,2000,5000,10000,20000,50000,100000,200000))

##  [1] (2e+04,5e+04] (1e+04,2e+04] (1e+04,2e+04] (1e+04,2e+04] (5e+03,1e+04]
##  [6] (5e+03,1e+04] (1e+04,2e+04] (1e+05,2e+05] (1e+04,2e+04] (2e+04,5e+04]
## [11] (2e+04,5e+04] (2e+03,5e+03] (1e+04,2e+04] (1e+04,2e+04] (2e+03,5e+03]
## [16] (5e+04,1e+05] (5e+03,1e+04]
## 7 Levels: (0,2e+03] (2e+03,5e+03] (5e+03,1e+04] ... (1e+05,2e+05]

cut 함수의 구분 기준(,)에 맞추어 heat.colors 함수를 연결하면 기준에 맞추어 R에서 자동으로 색상을 줍니다.

cut(confirmed, c(0,2000,5000,10000,20000,50000,100000,200000), heat.colors(7))

##  [1] #FFCC00 #FF9900 #FF9900 #FF9900 #FF6600 #FF6600 #FF9900 #FFFF80 #FF9900
## [10] #FFCC00 #FFCC00 #FF3300 #FF9900 #FF9900 #FF3300 #FFFF00 #FF6600
## Levels: #FF0000 #FF3300 #FF6600 #FF9900 #FFCC00 #FFFF00 #FFFF80

이제 kor1 자료에 색상을 넣은 후 이를 시각화 하겠습니다.

• 주의 : R이 색을 이해할 수 있도록 paste를 꼭 사용해야합니다.

사용하지 않으면 원하는 색으로 그림이 그려지지 않습니다.

kor1@data$col=cut(confirmed, c(0,2000,5000,10000,20000,50000,100000,200000), heat.colors(7))
plot(kor1, col=paste(kor1$col))

그림이 잘 그려졌지만 위험한 지역과 안전한 지역의 색이 반대로 된거 같습니다.

rev() 함수를 이용하여 색을 반대로 변환 하겠습니다.

kor1@data$col=cut(confirmed, c(0,2000,5000,10000,20000,50000,100000,200000), rev(heat.colors(7)))
plot(kor1, col=paste(kor1$col))

이제 그림의 중앙에 확인자 수를 넣어보겠습니다.

coordinates 함수를 이용하면 행정지역 별 가운데 위치를 R이 출력해줍니다.

coordinates(kor1)

##        [,1]     [,2]
## 1  129.0912 35.21142
## 10 127.8182 36.70061
## 11 126.8432 36.52545
## 12 128.5799 35.78117
## 13 127.3786 36.31868
## 14 128.2810 37.71682
## 15 126.8551 35.13442
## 16 127.1861 37.49094
## 17 128.7407 36.30965
## 2  128.2746 35.34559
## 3  126.6268 37.48770
## 4  126.5539 33.38431
## 5  127.1598 35.69758
## 6  126.9695 34.93609
## 7  127.2392 36.54758
## 8  126.9910 37.53330
## 9  129.2414 35.53645

확진자 수가 들어갈 위치의 x와 y의 위치를 lonlat으로 저장한 후 text 함수로 시각화 합니다.

legend(범례)를 추가하여 색상의 범위도 알려주겠습니다.

coordinates(kor1) -> lonlat
plot(kor1, col=paste(kor1$col))
text(x=lonlat[,1], y=lonlat[,2], labels=confirmed, cex=0.7)
legend("bottomright", 
       legend=c("0-2k","2k-5k","5k-10k","10k-20k","20k-50k","50k-100k","100k-200k"),
       fill=rev(heat.colors(7)), cex=0.6)

경기도와 서울특별시의 위치가 중복되기는 하지만 숫자가 들어갈 위치만 조정해주면 그림이 수정됩니다.

title을 이용하면 제목도 설정할 수 있습니다.

R을 이용하지 않더라도 ’SGIS 통계지리정보서비스’를 이용하면 지도에 시각화가 가능합니다.

SGIS 통계지리정보서비스

원하는 지역만 그리고 싶다면?

kor2는 시군구 그림입니다. 만약 수도권(경기도, 서울특별시, 인천광역시)만 선택하여 그림을 그리고 싶다면??

kor2@data를 통해 자료의 구조를 살펴보면 NAME_1을 이용하여 경기도, 서울특별시, 인천광역시를 쉽게 찾을 수 있습니다.

unique() 함수는 중복을 제외한 유일한 하나만 출력시켜주는 함수 입니다.

head(kor2@data)

##        GID_0      NAME_0   GID_1 NAME_1               NL_NAME_1     GID_2
## 184167   KOR South Korea KOR.1_1  Busan 부산광역시 | 釜山廣域市 KOR.1.1_1
## 184231   KOR South Korea KOR.1_1  Busan 부산광역시 | 釜山廣域市 KOR.1.2_1
## 184264   KOR South Korea KOR.1_1  Busan 부산광역시 | 釜山廣域市 KOR.1.3_1
## 184196   KOR South Korea KOR.1_1  Busan 부산광역시 | 釜山廣域市 KOR.1.4_1
## 184310   KOR South Korea KOR.1_1  Busan 부산광역시 | 釜山廣域市 KOR.1.5_1
## 184128   KOR South Korea KOR.1_1  Busan 부산광역시 | 釜山廣域市 KOR.1.6_1
##           NAME_2 VARNAME_2          NL_NAME_2 TYPE_2 ENGTYPE_2 CC_2 HASC_2
## 184167       Buk      <NA>         북구| 北區     Gu  District <NA>   <NA>
## 184231  Busanjin      <NA> 부산진구| 釜山鎭區     Gu  District <NA>   <NA>
## 184264      Dong      <NA>         동구| 東區     Gu  District <NA>   <NA>
## 184196   Dongnae      <NA>     동래구| 東萊區     Gu  District <NA>   <NA>
## 184310   Gangseo      <NA>     강서구| 江西區     Gu  District <NA>   <NA>
## 184128 Geumjeong      <NA>     금정구| 金井區     Gu  District <NA>   <NA>

unique(kor2@data$NAME_1)

##  [1] "Busan"             "Chungcheongbuk-do" "Chungcheongnam-do"
##  [4] "Daegu"             "Daejeon"           "Gangwon-do"       
##  [7] "Gwangju"           "Gyeonggi-do"       "Gyeongsangbuk-do" 
## [10] "Gyeongsangnam-do"  "Incheon"           "Jeju"             
## [13] "Jeollabuk-do"      "Jeollanam-do"      "Sejong"           
## [16] "Seoul"             "Ulsan"

이제 subset()를 이용하여 원하는 지역만 추출하겠습니다.

부산 아래의 섬 하나를 인천으로 인식하는 문제가 있습니다.

axis() 함수를 이용하여 그림의 범위를 살펴봅니다.

이 후 xlim과 ylim을 이용하여 그림의 범위를 지정하면 수도권 지역만 선책됩니다.

metro <- subset(kor2, NAME_1 %in% c("Seoul","Gyeonggi-do","Incheon"))
plot(metro)
axis(1)
axis(2)

plot(metro, xlim=c(125,128), ylim=c(36.5,38.5))

head(metro@data)

##        GID_0      NAME_0   GID_1      NAME_1       NL_NAME_1     GID_2
## 183025   KOR South Korea KOR.8_1 Gyeonggi-do 경기도 | 京畿道 KOR.8.1_1
## 183141   KOR South Korea KOR.8_1 Gyeonggi-do 경기도 | 京畿道 KOR.8.2_1
## 183010   KOR South Korea KOR.8_1 Gyeonggi-do 경기도 | 京畿道 KOR.8.3_1
## 182961   KOR South Korea KOR.8_1 Gyeonggi-do 경기도 | 京畿道 KOR.8.4_1
## 182758   KOR South Korea KOR.8_1 Gyeonggi-do 경기도 | 京畿道 KOR.8.5_1
## 182742   KOR South Korea KOR.8_1 Gyeonggi-do 경기도 | 京畿道 KOR.8.6_1
##             NAME_2 VARNAME_2          NL_NAME_2 TYPE_2 ENGTYPE_2 CC_2 HASC_2
## 183025       Ansan      <NA>     안산시| 安山市     Si      City <NA>   <NA>
## 183141     Ansoeng      <NA>     안성시| 安城市     Si      City <NA>   <NA>
## 183010      Anyang      <NA>     안양시| 安養市     Si      City <NA>   <NA>
## 182961     Bucheon      <NA>     부천시| 富川市     Si      City <NA>   <NA>
## 182758 Dongducheon      <NA> 동두천시| 東豆川市     Si      City <NA>   <NA>
## 182742    Gapyeong      <NA>     가평군| 加平郡    Gun    County <NA>   <NA>

위의 방법을 이용한다면 시군구 별 값과 색도 넣은 그림을 그릴 수 있습니다.

rgdal 패키지를 이용한 지도 시각화

GIS developer

위의 링크를 접속하면 TM 좌표계를 이용한 고용량의 시도, 시군구, 읍면동 자료를 다운로드 하실 수 있습니다.

shp 파일은 arcGIS에서 주로 사용합니다.

shp 파일은 rgdal 패키지의 readOGR 함수로 불러 올 수 있습니다.

shp를 처리에 관한 내용은 수업시간 강의를 진행하겠습니다.

전국표준노드링크도 readOGR 함수로 불러옵니다.

ggmap 패키지를 이용한 지도 시각화

ggmap을 이용하여 배경지도를 그린 후 결과를 시각화 할 수 있습니다.

if(!require(ggplot2))install.packages("ggplot2"); library(ggplot2)
if(!require(ggmap)) install.packages("ggmap"); library(ggmap)

ggmap(get_map(location='south korea', zoom=7))

을 실행시면 error가 발생합니다.

2018년까지 무료로 google지도를 R에서 다운로드 받아 시각화가 가능했지만,

지금은 API가 변경되어 카드를 등록하신 후 이용하실 수 있습니다.

구글맵 api key 발급방법

구글지도를 다운로드 받지 않더라도 아래 그림을 그릴 수 있습니다.

• 좌, 우, 위, 아래의 값을 지정해 주어야 지도 그림이 그려집니다.

library(ggmap)
kor <- c(left = 126.45, bottom = 33.49, right = 126.49, top = 33.52)
map <- get_stamenmap(kor, zoom = 14, maptype = "toner-lite")
ggmap(map)

map <- get_stamenmap(kor, zoom = 6, maptype = "watercolor")
ggmap(map)

만약 유럽 배경지도를 가지고 오고 싶다면?

europe <- c(left = -12, bottom = 35, right = 30, top = 63)
get_stamenmap(europe, zoom = 5) %>% ggmap()

한국지도를 배경으로 와이파이의 분포 그림을 그려보겠습니다.

kor <- c(left = 124, bottom = 33, right = 132, top = 39)
map <- get_stamenmap(kor, zoom = 7, maptype = "toner-lite")

wifi <- read.csv('F:/work/2020/1. 발표/R 시각화 (통계교육원)/day3/Chap13_map2/wifi.csv', header=T) #wifi.csv
head(wifi)

##   company      lat      lon
## 1      KT 37.74417 128.9056
## 2      KT 37.72806 128.9543
## 3      KT 37.75710 128.8900
## 4      KT 37.74769 128.8840
## 5      KT 37.74866 128.9073
## 6      KT 37.74281 128.8827

우리나라 이동통산시의 wifi 설치 위치들입니다.

geom_point를 이용하면 점으로 위치를 표현할 수 있지만 결과가 보기 좋지 않습니다.

ggmap(map) + geom_point(data=wifi, aes(x=lon, y=lat, color=company))

ggmap(map) + geom_point(data=wifi, aes(x=lon, y=lat, color=company)) + facet_grid(.~company)

stat_density_2d를 이용하면 2차원 밀도 그림(열섬그림)을 쉽게 그릴 수 있습니다.

ggmap(map) + stat_density_2d(data=wifi, aes(x=lon, y=lat))

ggmap(map) + 
  stat_density_2d(data=wifi, aes(x=lon, y=lat, fill=..level.., alpha=..level..), geom='polygon', size=2, bins=30)

..level..은 밀도의 수준을 계산하라는 함수이고 alpha는 투명도를 의미합니다.

geom=’polygon’에서 polygon은 다각형이라는 뜻입니다.

기본은 위에서 본 것처럼 선(線)으로 돼 있는데 그것 말고 도형으로 그리라고 명령을 준 겁니다.

선을 기준으로 size는 선 굵기, bins는 선 간격이라는 뜻입니다.

ggmap(map) + 
  stat_density_2d(data=wifi, aes(x=lon, y=lat, fill=..level.., alpha=..level..), geom='polygon', size=7, bins=28) +
  scale_fill_gradient(low='yellow', high='red', guide=F)+
  scale_alpha(range=c(0.02, 0.8), guide=F)

## Warning: It is deprecated to specify `guide = FALSE` to remove a guide. Please use `guide = "none"` instead.
## It is deprecated to specify `guide = FALSE` to remove a guide. Please use `guide = "none"` instead.

국내선 항공편도 그림을 그릴 수 있습니다.

airport <- read.csv('F:/work/2020/1. 발표/R 시각화 (통계교육원)/day3/Chap13_map2/airport.csv', header=T) #airport.csv
route <- read.csv('F:/work/2020/1. 발표/R 시각화 (통계교육원)/day3/Chap13_map2/route.csv', header=T) #route.csv

head(airport)

##   airport iata     lon     lat
## 1    강릉  KAG 128.944 37.7536
## 2    광주  KWJ 126.809 35.1264
## 3    군산  KUV 126.616 35.9038
## 4    김포  GMP 126.791 37.5583
## 5    대구  TAE 128.659 35.8941
## 6    목포  MPK 126.380 34.7589

head(route)

##   id airport     lon     lat
## 1  1     CJJ 127.499 36.7166
## 2  7     CJJ 127.499 36.7166
## 3 45     CJJ 127.499 36.7166
## 4 77     CJJ 127.499 36.7166
## 5  2     CJJ 127.499 36.7166
## 6  8     CJJ 127.499 36.7166

ggmap(map) + geom_point(data=airport, aes(x=lon, y=lat, size=3))

ggmap(map) + 
  geom_point(data=airport, aes(x=lon, y=lat, size=3))+
  geom_line(data=route, aes(x=lon, y=lat, group=id))

제주도 노선을 황금색으로 표현해 본다면?

CJU_route<-route[route$airport=="CJU","id"]
route1 <- route[route$id %in% CJU_route,]

ggmap(map) + 
  geom_point(data=airport, aes(x=lon, y=lat, size=3))+
  geom_line(data=route1, aes(x=lon, y=lat, group=id, size=2), col='gold')

leaflet 패키지를 이용한 대화형 그림

leaflat 패키지를 설치하면 오픈지도를 R에서 그릴 수 있습니다.

기본그림은 아래와 같습니다.

library(leaflet)
m <- leaflet() %>% 
  addTiles()
m

타일을 바꾸면 여러형태의 그림이 표현됩니다.

x <- leaflet() %>% 
   addTiles() %>% 
   setView( lng = 2.34, lat = 48.85, zoom = 5 ) %>% 
   addProviderTiles("NASAGIBS.ViirsEarthAtNight2012")
x

Nasa: NASAGIBS.ViirsEarthAtNight2012

Google map: Esri.WorldImagery

Gray: Esri.WorldGrayCanvas

Terrain: Esri.WorldTerrain

Topo Map: Esri.WorldTopoMap

더 알고 싶다면 here

m <- leaflet() %>% 
   addTiles() %>% 
   setView( lng = 2.34, lat = 48.85, zoom = 3 ) %>% 
   addProviderTiles("Esri.WorldImagery")
m

지도 자료도 시각화가 가능합니다.

# load example data (Fiji Earthquakes) + keep only 100 first lines
data(quakes)
quakes <-  head(quakes, 100)

# Create a color palette with handmade bins.
mybins <- seq(4, 6.5, by=0.5)
mypalette <- colorBin( palette="YlOrBr", domain=quakes$mag, na.color="transparent", bins=mybins)

# Prepare the text for the tooltip:
mytext <- paste(
   "Depth: ", quakes$depth, "<br/>", 
   "Stations: ", quakes$stations, "<br/>", 
   "Magnitude: ", quakes$mag, sep="") %>%
  lapply(htmltools::HTML)

# Final Map
m <- leaflet(quakes) %>% 
  addTiles()  %>% 
  setView( lat=-27, lng=170 , zoom=4) %>%
  addProviderTiles("Esri.WorldImagery") %>%
  addCircleMarkers(~long, ~lat, 
    fillColor = ~mypalette(mag), fillOpacity = 0.7, color="white", radius=8, stroke=FALSE,
    label = mytext,
    labelOptions = labelOptions( style = list("font-weight" = "normal", padding = "3px 8px"), textsize = "13px", direction = "auto")
  ) %>%
  addLegend( pal=mypalette, values=~mag, opacity=0.9, title = "Magnitude", position = "bottomright" )

m